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SPACE-TIME DOMAIN DECOMPOSITION AND MIXED 
FORMULATION FOR REDUCED FRACTURE MODELS 
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Abstract. In this paper we are interested in the ’’fast path” fracture and we aim to use global- 
in-time, nonoverlapping domain decomposition methods to model flow and transport problems in 
a porous medium containing such a fracture. We consider a reduced model in which the fracture 
is treated as an interface between the two subdomains. Two domain decomposition methods are 
considered: one uses the time-dependent SteklovPoincare operator and the other uses optimized 
Schwarz waveform relaxation (OSWR) based on Ventcell transmission conditions. For each method, 
a mixed formulation of an interface problem on the space-time interface is derived, and different 
time grids are employed to adapt to different time scales in the subdomains and in the fracture. 
Demonstrations of the well-posedness of the Ventcell subdomain problems is given for the mixed 
formulation. An analysis for the convergence factor of the OSWR algorithm is given in the case with 
fractures to compute the optimized parameters. Numerical results for two-dimensional problems 
with strong heterogeneities are presented to illustrate the performance of the two methods. 

Key words, mixed formulations, domain decomposition, reduced fracture model, optimized 
Schwarz waveform relaxation, Ventcell transmission conditions, time-dependent SteklovPoincare op¬ 
erator, convergence factor, nonconforming time grids 
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1. Introduction. In many simulations of time-dependent physical phenomena, 
the domain of calculation is a union of domains with different physical properties and 
in which the lengths of the domains and the time scales may be very different. In 
particular, this is the case for a domain where there exist fractures and faults. In 
such a case, the fluid flows rapidly through these paths while it moves much more 
slowly through the rock matrix. As a result, the contaminants present in the porous 
medium that travel with the fluid are transported faster than in the case when there is 
no fracture. Thus the time scales in the fractures and in the surrounding medium are 
very different, and in the context of simulation, one might want to use much smaller 
time steps in the fractures than in the rock matrix. For simplicity we consider the 
case in which the domain is separated into two matrix subdomains by a fracture. The 
permeability in the fracture can be larger or smaller than that in the surrounding 
medium. A large permeability fracture corresponds to a fast pathway and a small 
permeability fracture corresponds to a geological barrier. Here we are interested in the 
“fast path” fracture. Modeling flow in porous media with fractures is challenging and 
requires a multi-scale approach: first, the fractures represent strong heterogeneities 
as they have much higher or much lower permeability than that in the surrounding 
medium; second, the fracture width is much smaller than any reasonable parameter of 
spatial discretization. Thus, to tackle the problem, one might need to refine the mesh 
locally around the fractures. However, this is well-known to be very computationally 
costly and is not useful at the macroscopic scale (i.e. when the fractures can be 
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modeled individually). One possible approach is to treat the fractures as domains of 
co-dimension one, i.e. interfaces between subdomains (see ID 1313 M EZI ESI IH EH EZI 
and the references therein) so that one can avoid refining locally around the fractures. 
We point out that in these reduced fracture models, unlike in some discrete fracture 
models, interaction between the fractures and the surrounding porous medium is taken 
into account. 

We are concerned with algorithms for modeling flow and transport in porous 
media containing such fractures. In particular, in this article we investigate two 
space-time domain decomposition methods, well-suited to nonmatching time grids. 
We use mixed finite elements mm as they are mass conservative and they handle 
well heterogeneous and anisotropic diffusion tensors. 

The first method is a global-in-time preconditionned Schur method (GTP-Schur) 
which uses a SteklovPoincare-type operator. For stationary problems, this kind of 
method (see [101 0H j35] ) is known to be efficient for problems with strong hetero¬ 
geneity. It uses the so-called balancing domain decomposition (BDD) preconditioner 
introduced and analyzed in E3EZ1, and in m for mixed finite elements. It involves 
at each iteration the solution of local problems with Dirichlet and Neumann data and 
a coarse grid problem to propagate information globally and to ensure the consis¬ 
tency of the Neumann subdomain problems. An extension to the case of unsteady 
problems with the construction of the time-dependent Steklov-Poincare operator was 
introduced in |2H US ], where an interface problem on the space-time interfaces be¬ 
tween subdomains is derived. However, for the time-dependent Neumann-Neumann 
problems there are no difficulties concerning consistency, and we are dealing with only 
a small number of subdomains, so we consider only a Neumann-Neumann type pre¬ 
conditioner, an extension to the nonsteady case of the method of [33]. A Richardson 
iteration for the primal formulation was independently introduced in [231133] . and 
its convergence was analyzed. In the case of elliptic problems with fractures, a local 
preconditionner |2) significantly improves the convergence of the method. 

The second method is a global-in-time optimized Schwarz method (GTO-Schwarz) 
and uses the optimized Schwarz waveform relaxation (OSWR) approach. The OSWR 
and GTP-Schur methods are iterative methods that compute in the subdomains over 
the whole time interval, exchanging space-time boundary data through transmission 
conditions on the space-time interfaces. The OSWR algorithm uses more general 
(Robin or Ventcell) transmission operators in which coefficients can be optimized 
to improve convergence rates, see [2Qj 32[ [38]. The optimization of the Robin (or 
Ventcell) parameters was analyzed in [B] and the optimization method was extended 
to the case of discontinuous coefficients in 0 E ISl EH EH HE OH- Generalizations 
to heterogeneous problems with nonmatching time grids were introduced in Muni 
m m no nn 12 a mi m . More precisely, in mmm, a discontinuous Galerkin 
(DG) method for the time discretization of the OSWR algorithm was introduced and 
analyzed for the case of nonconforming time grids. A suitable time projection between 
subdomains is defined using an optimal projection algorithm as in mm with no 
additional grid. The classical Schwarz algorithm for stationary problems with mixed 
finite elements was analyzed in m ■ An OSWR method with Robin transmission 
conditions for a mixed formulation was proposed and analyzed in |28[ 129] , where a 
mixed form of an interface problem on the space-time interfaces between subdomains 
was derived. In m, an Optimized Schwarz method with Ventcell conditions in the 
context of mixed formulations was proposed. This method is not obtained in such a 
straightforward manner as in the case of primal formulations as Lagrange multipliers 
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have to be introduced on the interfaces to handle tangential derivatives involved in 
the Ventcell conditions. 

In this work, we define both a GTP-Schur and a GTO-Schwarz algorithm for a 
problem modeling flow of a single phase, compressible fluid in a porous medium with 
a fracture. A straightforward application of [25] would be to consider the fracture 
as a third subdomain and to take smaller time steps there. We consider instead 
however a reduced model in which the fracture is treated as an interface between two 
subdomains. 

The definition of the GTP-Schur method is a straightforward extension of that 
in [23- However, to define the GTO-Schwarz method, something more is needed: a 
linear combination between the pressure continuity equation and the fracture problem 
is used as a transmission condition (which leads naturally to Ventcell conditions), and 
a free parameter is used to accelerate the convergence rate. The well-posedness of the 
subdomain problems involved in the first approach was addressed in H3I221I35], using 
Galerkin’s method and suitable a priori estimates. In this paper, the proof of well- 
posedness of both the coupled model and the Ventcell subdomain problems involved 
in the GTO-Schwarz approach is shown to follow from a more general theorem that 
covers the two cases. 

Note that more general reduced models that can handle both large and small 
permeability fractures f55] introduce more complicated transmission conditions on the 
fracture-interface (in the form of Robin type conditions, where the Robin coefficient 
has a physical origin), and it is not yet clear how to formulate an associated domain 
decomposition problem with a parameter that can be optimized. 


This paper is organized as follows: in the remainder of the introduction (Subsec¬ 
tion we state an abstract existence and uniqueness theorem for evolution prob¬ 
lems in mixed form, the proof being deferred to Appendix |A} In Section[2j we consider 
a reduced model with a highly permeable fracture and prove its well-posedness. Then 
in Section[3]we consider the GTP-Schur approach, based on physical transmission con¬ 
ditions, for solving the resulting problem. Different preconditionners for this method 
are proposed. In Section [4] we consider the GTO-Schwarz method, based on more 
general (e.g. Ventcell) transmission conditions, for solving the resulting problem. We 
prove the well-posedness of the subdomain problems with Ventcell boundary condi¬ 
tions. In Section [5] we consider the semi-discrete problems in time using different 
time grids in the subdomains. Finally, in Section [fij results of two-dimensional (2D) 
numerical experiments comparing the different methods are discussed. 

1.1. Abstract evolution problems in mixed form. The goal of this section 
is to give an existence and uniqueness result for evolution problems posed in mixed 
form, in the spirit of the well-known theorem for weak parabolic problems (see for 
example PH vol. 5]). 

We consider two Hilbert spaces E, and M (.M will be identified with its dual), 
and assume we have continuous bilinear forms 

a : E x E —> R, b : E X M —> M, c : M x M —> R 

and a continuous linear form 


L : M 
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We study here an abstract version of a parabolic problem in mixed form: 

Find p G 7T 1 (0, T; M) and u G L 2 (0, T; E) such that, 

a(u,v) — b(v,p) = 0, \/v G E, 

(d t p,p) M + c(p,p) + b(u,p) = L(p), Vp, G M, (1.1) 

p(',0) = po, 

for some p 0 G M. 

We make the following hypotheses on the data: 

• The bilinear form a is positive definite on E: 

a(u,u) > 0 Vu G E, u / 0, (HI) 

so that a defines a norm on E, and we denote by E a the space E with the 
norm induced by the bilinear form a. Note however that this norm will not 
necessarily be equivalent to the initial norm on E. 

• The bilinear form c is positive semidefinite on M 

c(p,p) >0 Vp G M. (H2) 

• The bilinear forms a and b satisfy the following compatibility condition: there 
exists /? > 0 such that 

VugE, sup + ||u||| > ftlMlI- (H3) 

H&M WpWm 

• There exists a subspace W C M (with continuous embedding) on which the 
bilinear form b satisfies the stronger continuity property : there exists Ct > 0 
such that 


b(u,p) < C'b||u||E a ||M||w ) Vu G E and Vp G W. 


(H4) 


In most cases, the application of hypothesis 
form if it is written using the operator B : W 
form b, that is such that 


(H3) will appear in a more natural 
—> M associated with the bilinear 


\/u G W, p, G M, b(u , p) = ( Bu , p)m- 


Then, hypothesis (H3) can be written in the equivalent form: there exists fj > 0 such 
that 


Vu G E, ||Hu|| 


2 

m ' 


Ie 0 > 


U 


ll, 


(H3’ 


see the following remark. 

Remark 1.1. Hypothesis ( H31 is not equivalent to the inf-sup condition. The 
inf-sup condition expresses the surjectivity of B T , whereas here we need the ellipticity 
of B with respect to the norm defined by a. This also implies a form of compatibility 
between a and B. This implies that a is elliptic on the kernel of B, i.e. if Bu = 0, 
a{u,u) > /3||u|||. 


The basic existence and uniqueness result for problem (1.1) is the following: 


Theorem 1.2. Let M and E be Hilbert spaces, and let a, b and c be continuous, 
bilinear forms satisfying (HI) through (H3). Then , if L is a continuous linear form 
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on M and po £ W, where W C M satisfies (H4), then problem <o> has a unique 
solution, for which the following estimate holds: 


« 


L 2 (0,T;E) + 1 Ml °°(0,T;M) + 


'tP\\L 2 (0,T;M) < C(\\L\\ 


L 2 (0,T;M) 


Iboll 


( 1 . 2 ) 


The proof of the theorem will be given in Appendix [A] 

Remark 1.3. The case c = 0 is allowed, and is actually the most common case 


(cf Theorem 2.1). 


Remark 1.4. This result is a generalisation to the abstract setting of Lemma (3.1) 
in 1 35j). This problem has also been considered by Boffi and Gastaldi 1 '.12*) . but the es¬ 
timates given there (without proof) are different: they dispense with the regularity 
requirement po £ W, at the expense of introducing weighted estimates in time to cope 
with the possibility of a singularity at the initial time. A proof in a more general 
setting is given in W- However this proof uses semigroup theory (see Theorem f.l 
ofW> while the one we propose in this paper is with a priori estimates, in the same 
spirit as in M- 

We give a simple application of Theorem |1.2| (other applications will be given in 
Theorems 2.1 and |4.1| below). 


We consider the heat equation with Dirichlet boundary conditions in mixed form. 
For a domain C R d (d = 2 or 3) and T > 0, we look for p : LI x [0, T] —> R, 
solution of: 


dp A , 

m~ Ap = f 

p = 0 

p(x,0) =p 0 (x) 


in LI x [0, T] 

on dfl x [0, T] 
in Q. 


(1.3) 


To obtain the mixed form of ( |1.3| , we define the spaces E = H (div ,fi) and 
M = L 2 (Q), the bilinear forms a, b (here we will take c = 0) and the linear form L 


a : S x E 
b : E x M 


u ■ v 


R, a(u,v) = / 

Jn 

R, b(u,p)= / /rdivu, 
Jn 

L : M ^ R, L(p) = f fp. 

Jn' 


(1.4) 

(1.5) 

( 1 . 6 ) 


To apply Theorem 1.2 we check hypothesis (HI I to (H4) above. This is trivial 


for (HI) and (H2). To check (H3), we use the equivalent form (H3’|- Operator B is 


simply the divergence, so that 

\\Bu\\ 2 m + \\u\ || o = / I!divu|| 2 + / 
Jn Jn 


ur = u 


£) 


and (H3) is valid with (3 = 1. Last we check that (H4) is also valid with W = 
Using Green’s formula, we obtain 


b(u, p) = / pdivu = — u-Vp, 

Jn Jn 


from which (H4) follows. 
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2. A reduced fracture model. For a bounded domain f 1 of (d = 2, 3) with 
Lipschitz boundary <9fl and some fixed time T > 0, we consider the problem of flow 
of a single phase, compressible fluid written in mixed form as follows: 

sdtp + divu =q infix ( 0 ,T), 

u = —KVp infix (0,T), , 

p =0 on 9fl x (0, T), y ’ 

p(-, 0 ) =p 0 in fl, 

where p is the pressure, u the velocity, q the source term, s the storage coefficient and 
K a symmetric, time independent, hydraulic, conductivity tensor (see e.g. [28]). For 
simplicity we have imposed a homogeneous Diriclrlet condition on the boundary. 

We suppose that the fracture fl/ is a subdomain of fl , of thickness <5, that 
separates fl into two connected subdomains (see Figure [ 2 d] left where for visualization 
purposes the size of S is depicted as being relatively much larger than it is in reality), 

n\U f = fii ufi 2 , fiinfi 2 = 0- 

Also, for simplicity, we assume that fl/ consists of the intersection with fl of a line 
or plane 7 (depending on whether d = 2 or 3), together with the points x = x 1 + sn 
where x 1 € 7 , s £ (— |) and n is a unit vector normal to 7 . We denote by 7 * the 

part of the boundary of fli shared with the boundary of the fracture fl/: 

7 i = (<9fli n <9fl/) fl fl, * = 1 , 2 , 

and we denote by n, the unit, outward pointing, normal vector field on <9flj. We use 


,5 




Fig. 2.1. Left: The domain Q with the fracture Qf. Right: The domain Q with the interface- 
fracture 7 . 

the convention that for any scalar, vector or tensor valued function (j) defined on 
(pi denotes the restriction of (p to = 1,2,/. We rewrite problem 423 as the 
following transmission problem: 


sApi + div Ui 

= di 

in fli x ( 0 , T), 

f = 1,2, /, 

Ui 

= —KiS/pi 

in fli x ( 0 , T), 

i = 1,2, /, 

Pi 

= 0 

on ( 9 fli n < 9 fl) x (0, T), 

i = 1,2, /, 

Pi 

= Pf 

on 7 i x ( 0 , T), 

* = 1,2, 

u, ■ ni 

= Uf ■ Ui 

on 7 i x ( 0 , T), 

* = 1,2, 

Pi(-, o) 

= P0,i 

in fli, 

* = 1,2, /. 


In the reduced fracture model, the fracture fl/ is treated as a simple interface 
7 between subdomains fli and fl 2 (see Figure [271] right). We use the notation V r 
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(respectively div r ) for the tangential gradient (respectively tangential divergence) 
operators along the fracture 7 . We denote by s 7 and K 1 the storage coefficient and 
the permeability tensor in the (d — l)-dimensional fracture 7 . The reduced model 
that we consider was derived in [TJ 35]. It may be obtained by averaging across the 
transversal cross sections of the d-dimensional fracture fly. It consists of equations in 
the subdomains, 


sApi + div u, 

Ui 

Pi 

Pi 

Pit, 0 ) 


= Q i 

= - K.jS7pi 
= 0 
= Pi 
= Po,i 


in fij x (0, T), 
in fij x (0, T), 
on (dQi n <912) x (0,T), 
on 7 x (0,T), 
in f h, 


fori = 1,2, (2.3) 


and equations in the interface fracture 

s 7 <9 t p 7 + div T u 7 = («i • «!| 7 + u 2 • n 2 | 7 ) in7x(0,T), 

u~/ =—K 1 8'S7 T p 1 in7x(0,T), 

p 1 =0 on 97 x (0,T), 

Pit, 0) =Po,i in 7 . 


(2.4) 


These equations are the mass conservation equation and the Darcy equation in the 
subdomain together with the lower dimensional mass conservation and Darcy equa¬ 
tions in the fracture of co-dimension 1. These two systems are coupled: the fracture 
sees the subdomain through the source term in the conservation equation in the frac¬ 
ture which represents the difference between the fluid entering the fracture from one 
subdomain and that exiting through the other subdomain. Each subdomain sees the 
fracture through the Dirichlet boundary condition imposed on the part of its bound¬ 
ary common with the fracture. We make the hypothesis of the following compatibility 
conditions: p 0 ,i = Po ,7 011 7 , f° r * = 1,2. For a general mathematical treatment of 
this type of problem in the stationary case see m- 

To prove the well-posedness of problem (2.3)-(2.4), we shall use the abstract 


framework of Subsection 1.1 and apply Theorem 


1.2 We first write the weak for¬ 


mulation for problem (2.3)- (2.4), and define the appropriate function spaces, and 


the forms on these spaces. We use the convention that if V is a space of functions, 
then V is a space of vector functions having each component in V. For an arbitrary 
domain O , we denote by t, to the inner product in L 2 (0) or L 2 (0) and by and 
|| • ||o the L 2 {0)- norm or L 2 ((D)-norm. To write the weak formulation of ( |2.3| )-(2.4), 
we define the following Hilbert spaces: 

M = {p = (/ft,p- 2 ,/ft) e L 2 (VLi) x L 2 (fl 2 ) X L 2 ( 7 )} , 

E = {u = (v 1 ,v 2 l v 1 ) £ L 2 ( Hi) x L 2 (VL 2 ) x L 2 ( 7 ) : div v t £ L 2 (£li),i = 1,2, 


and div T u 7 — E Vi-n ih £ L 2 { 7 )}, 


equipped with the norms 


I 2 m 


= E IMn, + IKIIr 
2—1 

2 2 

IMII = E OHIfii + ll divt, illn 4 ) + IKII 7 + ll diy r w - E^-^dII 
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We define the following bilinear forms 


E x E 

(u,u) 

E x M 


((ir 7 <5) 1 u~ / ,v 


'7> t/ 7/ 


(u,n) 
c s : M x M —> 

(v,p) ^ 

and the linear form 


^ a(u,v) = Y ( K i 

i=1 

—s> K 

2 / 2 \ 

>->• b{u, n) ='Y (div Ui, Hi) u . + I div r u 7 - Y u * ' n d7> M7 ) > 


i =1 


i=l 


c s (? 7 ,m) = Y + (s 7 r? 7 ,/r 7 ) n 


L q : M 


^ Mm) = Y^'V^s 


With these spaces and forms, the weak form of (2.3)-(2.4) can be written as follows: 


Find p £ Ff 1 (0,T;M) and u £ A 2 (0,T;E) such that, 
a(u,v) — b(v } p) =0, Vue E, 
c s (d t p, m) + b{u, p) = L q (p), V/z £ M, 
together with the initial conditions 

Pi(-,0) =p 0 ,i in fli, 7 = 1,2, 

p 7 (-,0) =po, 7 in 7> 

for £ L 2 (flj), 7=1,2 and p 0 ,7 € A 2 (7)- We also define the space 

Hi := {^i = (Ati,/r 2 ,/r 7 ) € iT^Oi) x FV 1 (f2 2 ) x ^(7) : ft = 0 on n 3f2, 

and ft = ft on 7, 7 = 1,2}, 

equipped with the norm 

2 


(2.5) 

( 2 . 6 ) 


I Hi 


= Mm + Y 11^110, + IIV.M 7 II 


i=l 


The well-posedness of problem (2.5)-(2.6) is given by the following theorem: 

Theorem 2.1. Assume that there exist four positive constants s_ and S+, AT 
and K + such that 

• S- < Si(x) < s+ for a.e. x £ fli, 7=1,2, 

• s_ < s 7 (:r) < s + for a.e. x £ 7 , 

• <f T Kf 1 {x)q > A_|c| 2 , and \K z (x)<;\ < K + |c|, for a.e. x £ V? £ R d , 7 = 1,2, 

• q T (K 1 {x)5)~ 1 q > K_\r)\ 2 and |(A' 7 (:r)c)) _ 1 ? 7 | < K + \r]\ for a.e. x £ 7 , V 77 £ K d_1 . 
If q is in A 2 (0,T;M) and po = (Po,i,Mo, 2 ,Po, 7 ) in Hi then problem (2.5)-(2.6) has a 
unique solution ( p,u ) £ iT 1 (0,T;Af) x L 2 (0,T;E). 

Proof. First notice that under the assumptions on Sj and s 7 stated in the theorem, 
c s defines an inner product on M x M , and that the associated norm is equivalent to 
the original norm on M . 
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We will apply Theorem |1.2| in the case c = 0. The bilinear forms a and b are 
obviously continuous, and, with the hypotheses concerning Kj and AT 7 , a is positive 
definite on E. 


We now check hypothesis (H31. This is easiest to do using the operator form (H3’ 


where for this problem B is defined by 

/ 2 

Vu € E, Bu = I div Mi, div u 2 , div T u- — 


i=1 


The result follows from the hypothesis on A) and A' 7 and the definition of the norm 
on E. 

Last, to check hypothesis (H4), we can take W = H}, and use Green’s formula 
to see that for aeE and ^ £ W 

b{u , n) = ^2 (di v u i- A‘i)n i + div r u 7 - u t ■ n t | 7 , /z 7 

2—1 \ 2=1 / 7 
2 2 
= ^ ( [ u i , + (li; • flj| 7 , ^ — ( U 7 j ^tAG) 7 ~ ^ ' ( U 2 ’ W{| 7 , /* 7 ), 

2=1 

2 

= — ^ ' ( W 2 J ^b L i)o, i ~ ( u 7i 7 )^ 


2=1 


2=1 


because /ij| 7 = /x 7 . To conclude, we bound the terms of the right hand side: 
(Ui,V/Xj)Q. < A’+||tij||fl < 

(u 7 , V t /x 7 ) 7 < A+||u 7 || 7 || /x 7 ||_ffi( 7 ), 


from which hypothesis (H4| easily follows. □ 

It is natural to use domain decomposition methods for obtaining a numerical 
solution of problem ( 2 . 2 ) or problem (2.3)-(2.4), especially as these methods make 
it possible to take different time steps in the subdomains and in the fracture. For 
problem ( 2 . 2 ), it would be a straightforward application of the methods introduced 
in [29] while for problem (2.3)-(2.4), we need to derive a different formulation. In the 
following, we present two global-in-time domain decomposition methods for solving 
(2.3)-(2.41 based on different transmission conditions. A space-time interface problem, 
which will be solved iteratively, is derived for each approach. 

3. Global-in-time preconditioned Schur (GTP-Schur): using the time- 
dependent Steklov-Poincare operator. The Global-in-time preconditioned Schur 
(GTP-Schur) method is directly derived from the formulation of problem (2.3) - (2.4). 
To obtain the interface problem for this method, we need to introduce some notation. 
For a bounded domain O £ M. d (d = 2, 3) with Lipschitz boundary dO containing an 
open subset 7 C dO, we define the space 

Hl 7 (0) := {fi G H^O) : /* = 0 on (dO \ 7 )}. 

Then we define the following Dirichlet to Neumann operators 5j DtN , * = 1,2: 

<Sf tN : H\0,T;H^)) x A 2 (0, T; L 2 (^)) x L 2 ( 0 , T; (^( 7 ))' 

S? tN (A,g,p 0 ) l- t ^ ■ »j| 7 , 
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where ( Pi,Ui ), i = 1 , 2 , is the solution of the problem 


sApi + div Ut 

= Q 

Ui 

= -K 

Pi 

= 0 

Pi 

= A 

Pi{-, 0 ) 

= Po 


in flj x ( 0 , T), 
in fl> x ( 0 , T), 
on (<9fl; fl dfl) x (0,T), 
on 7 x ( 0 , T), 
in fl;. 


(3.1) 


Remark 3.1. A straightforward application of Theorem \1.S\ shows the well- 
posedness of subdomain problem & See also Ui3 i. | J5j for a direct proof. 

Problem (2.4) is reduced to an interface problem with unknowns A and rt 7 : 


Sjd t A + div r u 7 = 5f tN (A, g;,p 0 .i) in 7 x (0, T), 
w 7 = —RT 7 (5V r A in7x(0, T), 

A = 0 on c ?7 x (0, T), 

A(-, 0) =Po , 7 in 7 . 


(3.2) 


or equivalently 


s 7 d t A + div T Uj — X)i=i <5? tN (A, 0,0) = Ei=i ^ DtN (0i Qi,Po,i) in 7 x (0 ,T), 

u 7 =—K 1 SX7 r A in 7 x (0, T), 

A = 0 on dy x (0, T), 

A(-,0) =Po , 7 in 7, 

(3.3) 

or in compact form (space-time), 


5 


A 


= X- 


This problem is solved using an iterative solver such as GMRES since due to the time 
derivative the system is nonsymmetric. 

To improve the convergence of the iterative algorithm, we will consider two pre¬ 
conditioners. The first, introduced in [2], arises from the observation that the interface 
problem is dominated by the second order operator (div r (K^5\/ T )) since the Steklov- 
Poincare operator is of lower order (first order). This is even more the case when the 
permeability in the fracture is much larger than that in the surrounding domain. Thus 
one choice for a preconditioner is defined by taking the discrete counterpart of 
the operator (div T (i)f 7 (SV T )) _1 . We have 

Pf 0 c^ 2 (7) ~>£ 2 ( 7) 

£7 


where (p 7 ,it 7 ) is the solution of the problem 


div T u 1 = g 1 in 7 , 

u 1 = —K^S'Vrp^ in 7 , 
p-y = 0 on 97 . 

This preconditioner was introduced for elliptic problems, and it was shown numeri¬ 
cally [ 2 j that it significantly improves the convergence of the algorithm, especially, as 
mentionned before, for high permeability in the fracture. 
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A second possibility is to use the Neumann-Neumann preconditioner as was done 
in |28| 129) for ordinary domain decomposition algorithms (i.e. without fractures). 
The preconditioned problem is then 


with 


-p— 1 . ( - /oDtN\ —1 | /oDtN\ —1\ 

NN ' I. 0 l( < - > 1 ) +02(02 ) ), 

where ex.; : T x (0, T) —>• [0,1] is such that + tr 2 = 1. If K t = RJ and A; is constant 
in each subdomain then 

<Jj = -. 

Ai + S 2 

The operator (5P tN ) _1 , i = 1,2, is the inverse of the operator <Sf tN := Sf tN (-, 0,0), 
and is defined by 

(5f tN )“ 1 : L 2 (0, T; L 2 ( 7 )) -> H 1 (0, T; L 2 ( 7 )) 

( 4 DtN ) 1 (+) ^ Pi\~n 

where (p;,?*;), i = 1,2, is the solution of the problem 


sApi + div 

Ui 

= 0 

in fl; 

x (0, T), 


Ui 

= - K,,\7pi 

in 12 ; 

x (0, T), 


Pi 

= 0 

on (<9f2; n dQ) x ( 0 ,T), 

-u, ■ 

n, 

= <P 

on 7 

x (0,T), 

Pi(■, 

0 ) 

= 0 

in 



In Section | 6 j we will carry out numerical experiments and compare the performance 
of these two preconditioners. 


4. Global-in-time optimized Schwarz (GTO-Schwarz): using optimized 
Schwarz waveform relaxation. While the extension of the GTP-Schur method 
to handle the fracture model is straightforward, the extension of the GTO-Schwarz 
method to the fracture problem needs something more. Indeed, instead of impos¬ 
ing Dirichlet boundary conditions on 7 x (0, T) when solving the fracture problem 
as was done for the GTP-Schur method, for the GTO-Schwarz approach one uses 
optimized Robin transmission conditions. Thus, we introduce new transmission con¬ 
ditions, that combine the equation for continuity of the pressure across the fracture 
with the flow equations (2.4) in the fracture. These new transmission conditions con¬ 
tain a free parameter, which is used to accelerate the convergence. This is an extension 
of the OSWR method with optimized Robin parameters studied in [28] [29] in which 
Robin-to-Robin transmission conditions are considered in mixed form. Here how¬ 
ever, because of the fracture problem, we obtain what we will call Ventcell-to-Robin 
transmission conditions as described below. 


4.1. Ventcell-to-Robin transmission conditions. The new transmission con¬ 
ditions are derived by introducing Lagrange multipliers p, j7 , * = 1,2, with pj i7 repre¬ 
senting the trace on the interface 7 of the pressure pt in the subdomain Hj. As the 
pressure is continuous across the interface, one has 


Pi ,7 = P 2,7 = P 7 , on 7 x (0, T). 


(4.1) 
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We then rewrite the Darcy equation in the fracture associated with each pi ?7 as 

u 7j j = — Kj5\/ T pi yJ , on 7 x ( 0 , T), i = 1 , 2 . 

We have used the notation u 7 j j, i = 1,2, instead of u lr , to insist on the fact that u 7i j 
is not the tangential component of a trace of Ui on 7 . In fact, u 7 j j, i = 1,2, represents 
the tangential velocity in the fracture associated with the pressure Pi so that 


u 7) i = u 7 ,2 = « 7 , on 7 x ( 0 , T), * = 1 , 2 . 


With the notation introduced above, the flow equation (2.4) in the fracture can be 
rewritten, for * = 1,2, and j = (3 — i), as 


—Ut ■ rii + s 1 d t pi jl + div T u 7ii = —Uj rii, on 7 x (0, T), 

i* 7i j on 7 x (0 , X 1 ), 

Pi rf =0 on (I 7 x (0, T), 

Pvr(-,0) =Po, 7 1117. 


(4.2) 


In the context of domain decomposition, (4.1) and ( |4.2[ ) are the coupling conditions 
between the subdomains. As in the case without a fracture we take a linear com¬ 
bination of these conditions (for a parameter a > 0 ), but here we obtain equivalent 
Ventcell-to-Robin transmission conditions (instead of Robin-to-Robin): 


— «1 • Til + OPl j7 + SrydtPl ry + dfV T U-y^ 

u 7 ,l 

—u 2 ■ n 2 + ap 2, 7 + s 1 d t p 2 ,j + div T u 7i2 

W 7 ,2 

Using these transmission conditions together with the boundary and initial conditions 

Pi ,7 =P 2, 7 =0 on cby x (0, T), _ 

Pi, 7 (-, 0 ) =P 2 , 7 (-, 0 ) =P 0, 7 in 7 , 

the subdomain problem is obtained by imposing Ventcell boundary conditions on 
7 x (0, T), i = 1,2, j = 3 -i: 


-u 2 ■ n 1 + ap 2, 7 

T p\ n 


-ui • n 2 + api , 7 
-K 1 5V T p 2 , 1 


on 7 x (0,T), 

(4.3) 

0117 x ( 0 ,T), 


• n,: 


QPi,') 


SidtPi + div Ui 

= q 

in f l t x (0,T), 

Ui 

= -KiVpi 

in Qi x (0, T), 

0 j )7 + diV T Uryi 

— - Uj ' Tli CXPj^ry 

on 7 x (0, T), 


^f 7 ^VjPj j7 

in 7 x (0,T), 

Pi 

= 0 

on (9Uj n 9U) x (0,T), 

Pin 

= 0 

on c ?7 x (0,T), 

Pi{-, 0 ) 

= Po 

in fi-i, 

Pi,■> 0) 

0 

Oh 

II 

in 7 , 


(4-6) 

where the quantity on the right hand side of the third equation will be known in the 
context of an iterative method for solving (2.3)-(2.4). In the next subsection we prove 
that problem (4.6) is well-posed. 
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4.2. Well-posedness of the subdomain problem with Ventcell boundary 
conditions. For a bounded domain O C (d = 2, 3) with Lipschitz boundary dO 
containing an open subset 7 C 80, consider the following time-dependent problem 
written in mixed form with Dirichlet and Ventcell boundary conditions 


sodtpo + div uo 

= q 

in O x (0, T), 

u 0 

= —Ko^po 

in O x (0 ,T), 

uo ■ n + ap-y + Srydtp-y + div r u 7 

= e 1 

on dj x (0, T), 

lt 7 

= —KydS/rP-y 

in 7 x (0, T), 

Po 

= 0 

on (dO \ 7 ) x (0, T), 

V 7 

= 0 

on 97 x (0, T), 

Po( m , 0) 

= Po,o 

in fli, 

p 7 (-, 0 ) 

II 

O 

in 7 , 


where 0 7 is a function defined on 7 x (0, T), and a £ R, a > 0. In order to write the 
weak formulation of ( |4.7| ), we need to define the following Hilbert spaces: 

M 0 = {p = (p a ,P-y) G L\0) x L 2 ( 7 )} , 

S 0 = {t> = (vo,v 7 ) £ L 2 (0) x L 2 ("f) : divi>o G if(div ,0) and (div T Vy —vq ■ «| 7 ) G L 2 
equipped with the norms 

Mmo = \Wo\\1 + \\^ 

Hilo = l|uo||o + ||divVo||o + IKII 7 + ||div r t> 7 - t> 0 ■ n, 7 1| 2 . 


Then define the bilinear forms 


ao : 

Do x So 

—y R, 

ao(u,v ) = 

(K 0 1 uo,vo^ ) 

1 O (C^ 7 ^) ^77^7) 

bo : 

So x Mq 

—^ R, 

bo(u, p) = 

(div Uo,po)o + (div T Uy —u 0 ■ n\ 

co : 

Mq X Mq 

—^ R, 

c o(Vi m) = 

(cc?7 7 , /r 7 )^ , 


c s ,o ■ 

Mq X Mq 

— y R, 

c s ,o(v,p) = 

(sovo, Po)o 

T P^) y > 

and the linear form 






Lq ,0 '■ 

Mq — 

g L qt o(p) 

= (q, vo)o + 

(9y, Py)y ■ 


With these spaces and forms, the weak form of ( |4.7| can be written as follows: 
For a.e. t £ (0 ,T), find p(t) £ Mq and u(t) £ So such that 


ao(u,v) - bo(v,p) =0 

Vv £ So, 

(4.8) 

c s ,o (d t p, p) + co (p, p) + b a (u,p) = L q>0 (p) 

V/i £ M 0 , 

together with the initial conditions 



Po(-, 0) =Po,o in O, 
Py(-,0) =Po , 7 in 7- 


(4.9) 


We will also make use of 


7 ) ~ Hl„(P) x Hl{ 7 ). 


The well-posedness of problem (4.8 1-(|4.9[) is given by the following theorem: 


Theorem 4.1. Assume that there exist positive constants s_, S+, K_, K+ with 
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• S- < so{x) < s+ for a.e. x G O, 

• s_ < s 7 (x) < s + for a.e. x G 7 , 

• q T K~^(x)<; > AV|<r| 2 , and \Ki(x)<;\ < K + \<;\, for a.e. x G O, V? G R d , 

• if 1 K~ 1 \x)5p > A'_|ry| 2 and \(K 1 (x)S)~ 1 ri\ < K + for a.e. x G 7 , V77 G 

If q is in A 2 (0, T; M@), p 0 in ,7) and 0 7 in L 2 (0,T; L 2 ( 7)) t/ien problem (4.8)- 

(4.9) /las a unique solution (p,u) G H l (f),T\Mo) x A 2 (0,T;£e>). 


d— 1 


Proof. As in the case with Dirichlet boundary conditions, we apply Theorem 1.2 


Again notice that, under the hypotheses on so and s 7 , c Si e> defines an inner product 
on M, equivalent to its usual inner product. 

• It is clear that a, b and c are all continuous forms, that a is positive definite 
(due to the hypotheses on Kq and AT 7 ) and that cq is positive (since a > 0), so that 


hypotheses (HI) and (H2) in Theorem 1.2 hold. 


To verify hypothesis (|H3|), we define the operator B by 


Bu = (div uo, div r m 7 — uo ■ n 


l \i) 


NIL 


and it follows easily that 

||-B«||m + ao(u,u) = ||div« 0 ||o + ||div T u 7 -u 0 -n | 7 || 7 

+ (Kq 1 u 0 ,u 0 ) 0 + ({K 1 5)~ 1 u 11 u 1 ) i > 

for some (3 > 0, again because of the lower bounds on Kq and AT 7 . 

• Last, to check (H4), we proceed as in Theorem 2.1 and use Green’s formula for 
uGS and p G 

b 0 {u,n) = (div Uo,Ho)o + ( div r ^7 ~ u o •n| 7 ,/z 7 ) 7 

= — ( UQ) V+ (uq '^ 17 ? — (^7? V — (uq • W|^y, M7)^ 5 

from which the proof of the theorem follows. □ 

4.3. The interface problem. As for the GTP-Schur method, we derive an 
interface problem which in this case is associated with Ventcell-to-Robin transmission 
conditions (4.3)-(4.4). Towards this end, we define the following Ventcell-to-Robin 
operator <S^ tK , which depends on the parameter a, for i = 1, 2; j = (3 — i ): 

^ vtR : L 2 (0,T; A 2 ( 7 )) x A 2 (0, T; A 2 (^)) x x -+ A 2 (0, T; L 2 ( 7 )) 

5, vtR (6» 7 , Q,Po,Po,j) '-t -«i • nj| 7 + ap i)7 , 

where (pi,u 2 ,pi rf .u~ ul ) is the solution of the subdomain problem with Ventcell bound¬ 
ary conditions 


sApi + div Ui 

= q 

in fl,; x (0, T), 

Ui 

= -K,Vp t 

in H,; x (0, T), 

Pi , 7 + div r w 7:i 

= 0 7 

on 7 x ( 0 ,T), 

^7 ,i 

— A" 7 ^V rPij'Y 

in 7 x (0, T), 

Pi 

= 0 

on (9Hi n 9H) x (0, T), 

Pi ,7 

= 0 

on 97 x (0,T), 


= Po 

in Hi, 

Pi, 7 O» d ) 

il 

O 

in 7 . 


(4.10) 
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The interface problem with two Lagrange multipliers is then 

0 7 ,i = 5j tR (0 7j2 ,92,Po, 2 ,^ 0 , 7 ) , 

,, iD , on 7 x ( 0 , 

0 7 ,2 =5^ tR (0 7) i,gi,po,i)Po, 7 ) 


(4.11) 


or equivalently 


0 7il -5 2 v «(0 x2 ,0,0,0) 

0 7 ,2-Sy tR (0 7l i, 0,0,0) 


= 5y tR (0,g 2 ,Po,2,Po, 7 ) 
= ^y tR (0, qi , p 0 ,i; Po, 7 ) 


on 7 x (0, T). (4.12) 


The discrete counterpart of this problem can be solved iteratively using Jacobi iter¬ 
ations or GMRES. The former choice yields an algorithm equivalent to the OSWR 


algorithm for the reduced fracture model (2.3) - (2.4) and is written as follows: start¬ 


ing with an initial guess 0° •, j = 3 — *, on 7 x (0, T) for the first iteration, 


-u) 


■m 


a p. 


2,7 


y %>° 7 + div T 


u =6^ 
“ 7> * 


then at the k th iteration, k = 1 ,..., solve in each subdomain the time-dependent 
problem, for * = 1,2; j = (3 — i), 




■ n, 


aPi, 7 


sidtPi 

Apt 


cliv. 


div u k 

u\ 

u i,i 
U j,i 

Pi 
Pi , 7 

pH-, 0 ) 

k 


J? 

2,7 


in f2j x (0, T), 
in 0* x (0, T), 
on 7 x (0, T), 
on 7 x (0, T), 
on (<9fl; n dfl) x ( 0 ,T), 
on dj x (0, T), 
in fli, 
in 7 , 

(4.13) 

k ~ 1 on 7 x (0, T). 

The convergence of algorithm ( |4.13 ) depends on the choice of the parameter a. 
Thus we extend the analysis for the convergence factor of the OSWR algorithm derived 
in the case without fractures iinama to this algorithm and from that, one can 
calculate the optimal or optimized values of the parameter a. 


Pi, A A 


= <H 

- —KiWPi 

_ fik— 1 

7 ,3 

= -K f , T 8V T p: 

= 0 
= 0 

= P0,i 

= ^0,7 


with e k ~H = -uH 1 


n 


+ ap JtJ 


4.4. Convergence factor formula for computing the optimized parame¬ 
ter. In this section, we extend the two domain analysis HHSlUiGIllMI to derive the 
convergence factor of the OSWR algorithm introduced in Section |4.3| for a reduced 
fracture model for compressible flow. Towards this end, we consider the two half¬ 
space decomposition = R _ x R, £l + = R + x R and write the OSWR algorithm, 
applied to the fractured model, in the primal formulation: at the k th Jacobi iteration, 
solve 


K 


dp^_ 

d 


s-d t p k _ + div {-KS7p k _) 
+ apt + s-ydt.p k _ + div T (— Kf^SS/rpt) 


P 1 K-,0) 


= q 


= K 


dp 


k -1 


dri- 


= Po 


in S2_ x (0, T), 
+ap k + ~ 1 

on 7 x (0, T), 
in f 2 _, 


(4.14) 
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and 


K 


dp\ 

dn+ 


s+d t p% + div (-iWVp+) = q 

+ ap\ + Sl d t p\ + div r (-Kf,rSV T p+) = K. 


dpt 


k -1 


in fi_|_ x (0, T), 


dn^ 


■+ap k _ 1 


on 7 x ( 0 , T), 

P+0,0) =Po in n+, 

(4.15) 

where j = {x = 0} is the fracture. We assume that the permeability is isotropic: 

Kpm — A±I, and K — Af : 

where / is the 2D identity matrix, and that the solution of the problem decays at 
infinity. As the problem is linear, we only consider q = 0 and po = 0, and analyse 


the convergence of (4.14)-(4.15) to the zero solution. We use a Fourier transform 


in time and in the y direction with parameters u> and y, respectively, to obtain the 
Fourier functions p± in time t and y of p± , as the solutions to the ordinary differential 
equation in x 


+ ( siu} +^ 2 ) p = °- 


Thus 


p = A(y,u J )e r+x + B(y,u J )e r ~t 
where 7 -± are the roots of the characteristic equation 

— Ar 2 + (siw + Ay 2 ) = 0, 


so 


? ,=t = ±-, A = 4A (siuj + Ay 2 ) . 

2A v ’ 

Here and throughout this article, we use the square root symbol y to denote the 
complex square root with positive real part. In order to work with at least square 
integrable functions in time and space, we look for solutions which do not increase 
exponentially in x. Since 3?? ,+ > 0 and 5ir - < 0, we obtain 

p k _ = A k (y,u])e r+( - s ~’ si -’ r >’ u '> x , 

p\ = B k {y,uj)e r ~( s +’ Si +’ r i^) x . 


Substituting these formulas into the transmission conditions on the interface 7 x (0, T) 
(i.e. the second equations of (4.141 and (4.15)), we find 


(A_?’ + (s_, A-,y,uj) + a + s^iuj + AfSy 2 ) p^(0,y,u}) 

= (A + r~{s + ,A +1 y,uj)+a)p k + - 1 (0,y,uj), 

(—A + r~(s + ,A + ,y,u}) + a + s^ioj + A f 5y 2 ) p\{ 0,y,u) 

= (-J?_r+(s_,il_,77,w) + a)p* _1 (0,77 ,w). 

(4.16) 
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From (4.16) using induction ans using £ to denote s 7 *w + M.fSr/ 2 , we obtain 
&+r~ (s + , A + , 77 , w) + a 


pL k (0,v,u) = 


A_r+(s_,£_,*7,w) + a + ( 
M + r~(s + ,& + ,r/,uj ) + a 
,&-r+(s-, &-,T),u.>) + a + £ 
= PfP-(0,T},w), 


K (0,rj,v) 


—A_r + (s_, A_, p, oj) + a 
-8. + r-(s+,M. + 1 r),uj) + a + C 


f_ k ~ 


where 


/ R + r (s+,&+,r],uj) + a \ f — A_r + (s_,A_, 77,1c) + a \ 
V^-r+(s_,^_,77,w) + a + Cj \~&+r~ (s+, &+,r],u;) + a + ()) ’ 


is the convergence factor of the algorithm ( |4.14| )-( |4~T5 ). Similarly, we obtain 

p 2 + (0,r),u) = p k f p° + ( 0,r],uj), 


Thus, we can calculate the parameter a in such a way as to minimize this continuous 
convergence factor: 


min 

a>0 



max |pf(s + ,& + ,s_,£_,a:,f7,w)| , 


(4.17) 


where L is the length of the fracture, h is the spatial mesh size, T is the final time 
and At is the maximum time step of the discretization in time. 

Remark 4.2. One could make use of the two-sided Robin as in In this 

article, the optimized one-sided Robin parameter works well since in the test case we 
considered, the two subdomains ffi and O 2 (representing the rock matrix) have similar 
physical properties (though a comparison of the performance of the one-sided and two- 
sided Robin might be considered). In our applications, the fracture is assumed to 
have much larger permeability than the surrounding domain, which suggests that the 
time step inside the fracture should be small compared with that of the surrounding 
matrix subdomains. As both of the methods derived in Sections [3] and [4] are global 
in time, i.e. the subdomain problem is solved over the whole time interval before 
the information is exchanged on the space-time interface, we can use different time 
steps in the fracture and in the rock matrix. In the next section, we consider the 
semi-discrete problem in time with nonconforming time grids. 


5. Nonconforming discretization in time. Let 71,72 and 7^ be three dif¬ 
ferent partitions of the time interval (0,T) into subintervals J l m = {t l m _ 1 ,t l rn ) for 
to = 1,..., Mj , and * = 1,2, 7 , (see Figure 5.1). For simplicity, we consider uniform 
partitions only, and denote by At*, * = 1 , 2,7 th e corresponding time steps. Assume 
that A t 1 <C A ti, i = 1,2. We use the lowest order discontinuous Galerkin method 
m [251 145 ) , which is a modified backward Euler method. The same idea can be 
generalized to higher order methods. 

We denote by To (71, L 2 ( 7 )) the space of functions piecewise constant in time on 
grid 71 with values in L 2 ( 7 ): 


Po(71,L 2 (7)) = {if : (0, T) -A- L 2 ( 7 ),-0 is constant on J, VJ S 71} . 

In order to exchange data on the space-time interface between different time grids, 
we use, for i,j in { 1 , 2 , 7 }, the L 2 projection Iljj from P 0 (71,L 2 ( 7)) to Po{Tj , L 2 ^)): 
for if G Pq(71, P 2 (7)), Wjiif\jj is the average value of if on J^, for to = 1,..., Mj. 


2 (0 ,r],u) 
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A U 


T = M 1 At 1 = M 2 At 2 = M 7 At-y 


Fig. 5.1. Nonconforming time grids in the rock matrix and in the fracture. 


For the GTP-Schur method. The unknown A on the interface represents 
the fracture pressure, thus A is piecewise constant in time on grid Py. In order to 


obtain Dirichlet boundary data for solving subdomain problem (3.1), we project A 
into P 0 (Ti,L 2 ( 7 )), for * = 1,2: 

Pi = n i 7 (A), on 7 , * = 1,2. 


The semi-discrete counterpart of the interface problem (3.2) is obtained by weakly 
enforcing the fracture problem over each time sub-interval of Py as follows 

rp +1 r t ™ +1 / 2 \ 

s 7 (A m+1 - A m ) + div T < +1 = («S° tN (n^(A), %,p o ,0) J, 


u™ +1 = -Ky6X7 T X m+1 

A m+1 = 0 


A 0 = Po ,7 


m 7, 

ondy, (5.1) 
in 7 , 


where A m = A , for m = 0,..., M 1 — 1. 

For a function piecewise constant in time p on the fine grid Py, the semi-discrete 
Neumann-Neumann preconditioner (still denoted by P~^ N ) is defined by: 


Pnn<p ■■= E ( (s? tN ) 1 ( n *7 


(5.2) 


where we have solved the subdomain problem with Neumann-Neumann data projected 
from Py onto 71, i = 1 , 2 , then extracted the pressure trace on the interface and 
projected backward from 71 onto Py. Thus the interface problem is defined on the 
fracture time grid. 


Remark 5.1. In the nonconforming semi-discrete (in time) case, we see from (5.2) 
that Pf/jn is not strictly a preconditioner and may affect the accuracy of the scheme 
due the projection operators used to define it. This is indeed what we observed in the 
numerical experiments in Section [d| 

For the GTO-Schwarz method. In the GTO-Schwarz method, there are two 
interface unknowns representing the linear combination of the fracture pressure and 
some terms from the fracture problem. Thus we let 0 7i j £ Po(Py, L 2 (T)), for j = 1, 2. 


In order to obtain Ventcell boundary data for solving the subdomain problem (4.7) 
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we project 0 7J onto the %, for i = 1, 2; j = 3 — i : 

—Ui ■ n, + api + s-ydtPi + div r u 7 = IIj 7 (0 7i j), on 7 , i = 1, 2. 


Remark 5.2. This setting is different from the case of usual domain decom¬ 
position (without fractures) analyzed in (281 (29 j/, where the two interface unknowns 
represent the Robin data in each subdomain and thus are chosen to be constant on the 
associated subdomain’s time grid, i.e. 6~i £ PoOl, L 2 (T)), fori = 1,2. 


The semidiscrete-in-time counterpart of (4.11) is weakly enforced over each time 


sub-interval of the fracture time grid as follows: for all m = 0,..., M 1 — 1, 

J # 7,1 =J n 7 2 (Sy iR (n 2l (0 1 ,2),q2 i PO,l,PO,~ / )) , 


7 7 ,2 


n 7 i ( 5 j /tR (IIi 7 ( 0 7 7 ), qi,Po,i,Po,j)) ) 


0117 , 


(5.3) 


Remark 5.3. We point out that with the GTO-Schwarz method as with the 
GTP-Schur method preconditioned by a Neumann-Neumann preconditioner (cf. Re¬ 


mark 5.1), we can not hope to gain in accuracy in the fracture by using a finer grid 


in the fracture only since the fracture problem is actually solved on the coarser time 
grids of the two subdomains. We will see this in the numerical experiments. 

6. Numerical results. In all of the numerical experiments, for the spatial dis¬ 
cretization we use mixed finite elements with the lowest order Raviart-Thomas spaces 
on rectangles Pi. 

Remark 6.1. The subdomain problem of the GTO-Schwarz method corresponding 
to Ventcell boundary conditions is somewhat more complicated than that of GTP- 


Schur method (problem (3.1 )). Consequently, for solving problem (4.6 1 , one needs 


to introduce Lagrange multipliers (see e.g. ' flJ[ \ 44V on the interface to handle the 
Ventcell conditions (representing the fracture problem). 

We carry out some preliminary experiments to investigate the numerical perfor¬ 
mance of the two methods proposed above. We consider the test case pictured in 
Figure |6.1| where the domain is a rectangle of dimension 2x1 and is divided into 
two equally sized subdomains by a fracture of width S = 0.001 parallel to the y 
axis. The permeability tensors in the subdomains and in the fracture are isotropic: 
K = , i = 1,2, /, and Jl, is assumed to be constant. Here we choose = £2 = 1 

and M.f = 10 3 (so that M.f5 =1). A pressure drop of 1 from the bottom to the top 
of the fracture is imposed. On the external boundaries of the subdomains a no flow 
boundary condition is imposed except on the lower fifth (length 0 . 2 ) of both lateral 
boundaries where a Dirichlet condition is imposed: p = 1 on the right and p = 0 on 


the left. See Figure 6.1 


We consider a uniform rectangular mesh with size h = 1/100. In time, we fix 
T = 0.5 and use uniform time partitions in the subdomains with time step At,, i = 1,2, 
and in the fracture with varying time step At 7 . We first consider the case with the 
same time step throughout the domain, Ati = A t 2 = A t 1 = A t = T/ 300. 

In Figure [672] snapshots at different times of the pressure field and of the flow field 
(on a coarse grid for visualization) are shown. The length of each arrow is proportional 
to the magnitude of the velocity it represents and the red arrows represent the flow 
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Fig. 6.1. Geometry of the test case where the fracture is considered as an interface. 


in the fracture. We see that the flow field is a combination of flow in the fracture and 
flow going from right to left in the rest of the porous medium and there is interaction 
between them as some fluid flows out of the fracture (near the bottom) and some 
flows into it (near the top at later times). Since Sif = 1,2, the velocity is 

much larger in the fracture than in the surrounding medium. 

Next, in order to analyze the convergence behavior of both methods, we con¬ 
sider the problem with homogeneous Dirichlet boundary conditions (i.e., the solu¬ 
tion converges to zero). We start with a random initial guess on the space-time 
interface-fracture and use GMRES as an iterative solver and compute the error in the 
L 2 (0, T; L 2 (S2))-norm for the pressure p and for the velocity u. We stop the iteration 
when the relative error is less than 10 -6 . We consider four algorithms: the GTP-Schur 
method with no preconditioner, the GTP-Schur method with the local preconditioner, 
the GTP-Schur method with the Neumann-Neumann preconditioner and the GTO- 
Schwarz method with the optimized Robin parameter. We compare the convergence 
behavior of these four algorithms in terms of the number of iterations. Note however 
that for the GTP-Schur method with the Neumann-Neumann preconditioner the cost 
per iteration is roughly twice as large as that of the other methods. 


In Figure [6731 the error curves versus the number of iterations are shown: the error 
in p (on the left) and in u (on the right). We see that the GTP-Schur method with no 
preconditioner (the blue curves) converges extremely slowly (after 500 iterations, the 
error, both in p and in u, is about 10 _1 ). The performance of the GTP-Schur method 
with the local preconditioner (the green curves) is better but still quite slow: it requires 
about 350 iterations to reach an error reduction of 10 -6 . The Neumann-Neumann 
preconditioner (the cyan curves) further improves the convergence rate about 150 
iterations are needed to obtain a similar error reduction. The GTO-Schwarz method 
needs only 6 iterations to reduce the error to 10~ 6 and thus the convergence of the 
GTO-Schwarz method is much faster than that of the other algorithms (at least by a 
factor of 25). This is due to the use of the optimized parameter a. In Figure 6.4 we 


show the error in u (in logarithmic scale) after 10 Jacobi iterations for various values 
of a. We see that the optimized Robin parameter (the red star) is located close to 
those giving the smallest error after the same number of iterations. Also we observe 
that the convergence can be significantly slower if a is not chosen well. 

Next, we study the behavior of three of the algorithms when nonconforming 
time grids are used. For this we again use the nonhomogeneous boundary conditions 
depicted in Figure [O] In all cases, we consider equal time steps for the subdomains 
as they have the same permeability: Ati = At 2 = At m . We examine three time grids 
as follows: 


• Time grid 1 (conforming coarse): At m = Atf = T/ 100. 

• Time grid 2 (nonconforming): At m = T /100 and At/ = T/500. 
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Fig. 6.2. Snapshots of the pressure field (left) and flow field (right) at t = T/ 300, t = T/4, 
t = T/2 and t = T respectively (from top to bottom). 
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Fig. 6.3. Convergence curves for the compressible flow: errors in p (on the left) and in u (on 
the right) - GTP-Schur method with no preconditioner (blue), with local preconditioner (green) and 
with Neumann-Neumann preconditioner (cyan) and GTO-Schwarz method (red). 



Fig. 6.4. L 2 velocity error (in logarithmic scale) after 10 Jacobi iterations for various values 
of the Robin parameter. The red star shows the optimized parameters computed by numerically 
minimizing the continuous convergence factor. 


• Time grid 3 (conforming fine): A t m = Atf = T/ 500. 

We start with a zero initial guess on the space-time interface and stop the GMRES 
iterations when the relative residual is less than 10 -6 . In Figure 6.5 we show the 
relative residual versus the number of iterations for three schemes: the GTP-Schur 
method with the local preconditioner, the GTP-Schur method with the Neumann- 
Neumann preconditioner and the GTO-Schwarz method with an optimized Robin 
parameter. We see that the GTO-Schwarz method still performs better than the 
GTP-Schur method, and the GTP-Schur method with the Neumann-Neumann pre¬ 
conditioner still converges faster than with the local preconditioner. The convergence 
rate of both the GTP-Schur method with the Neumann-Neumann preconditioner and 
the GTO-Schwarz method are almost independent of the time grid (the number of 
iterations does not change with the time grid) while that of the local preconditioner 
significantly depends on the temporal grid in the fracture. We also notice that the 
behavior of all three methods in the cases of nonconforming and conforming fine grids 
are very similar. 
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Fig. 6.5. Relative residual with GMRES for different time grids: GTP-Schur method with the 
local preconditioner (green), GTP-Schur method with the Neumann-Neumann preconditioner (cyan) 
and GTO-Schwarz method (red). 


Now we analyze the error (in time) of the three algorithms for each of the three 
time grids. A reference solution is obtained by solving problem (2.3) - (2.4) directly 
on a very fine time grid At = T /2000. The L 2 — L 2 error of the difference between 
the multi-domain and the reference solutions at each iteration is computed. We 
distinguish two different errors: error in the rock matrix L 2 (0, T; L 2 (flj)), * = 1,2, 
and error in the fracture L 2 (0, T; L 2 fy)). Figures 6.6 and 6.7 show the pressure error 
in the subdomains and in the fracture respectively. 


We first observe that the error in the subdomains after convergence (Figure 6.6) 
in the nonconforming case (Time grid 2) is equal to that in the conforming coarse 
case (Time grid 1) for all three algorithms. This is as expected as we use the same 
time step At m = T /100 in the matrix for both of these grids. However, as already 
pointed out in Remark |5.3[ though one might hope that the error in the fracture 
(Figure |6.7[ ) in the nonconforming case is close to that in the conforming fine grid 
case (Time grid 3), this can only be the case for GTP-Schur method with the local 
preconditioner. Only for this case do we actually solve the fracture problem on the 
fine grid. For the other algorithms, the fracture error of the nonconforming case is 
equal to that of the conforming coarse grid instead (see Remark 5.11. However, none 
of the methods deteriorates the accuracy because of nonconforming time grids. 

Remark 6.2. While the GTO-Schwarz method does not make it particularly 
useful to use a finer time grid in the fracture, it does give a rather remarkable con¬ 
vergence speed. For the advection-diffusion problem with an explicit time scheme for 
advection, one of the main advantages of using smaller time steps in the fracture is to 
avoid imposing a time step in the two subdomains dictated by the CFL number of the 
equation in the fracture. Thus we are hopeful that this algorithm will be useful when 
coupled with the advection equation simply for the convergence speed that it gives. We 
add however that we are still pursuing some ideas for modifying this scheme to obtain 
an algorithm that can take advantage of smaller time steps in the fracture for the 
diffusion equation. 


Conclusion. We consider two domain decomposition methods for modeling the 
compressible flow in fractured porous media in which the fractures are assumed to 
be much more permeable than the surrounding medium. Two space-time interface 
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GTP-Schur - local precond. GTP-Schur - NN precond. GTO-Schwarz 


Fig. 6.6. L 2 pressure error in the rock matrix: Time grid 1 (blue), Time grid 2 (magenta), 
Time grid 3 (black). 




GTP-Schur - local precond. GTP-Schur - NN precond. 


GTO-Schwarz 


Fig. 6.7. L 2 pressure error in the fracture: Time grid 1 (blue), Time grid 2 (magenta), Time 
grid 3 (black). 


problems are formulated using the time-dependent Dirichlet-to-Neumann and the 
Ventcell-to-Robin operators respectively, so that different time discretizations in the 
subdomains and in the fracture can be adapted. For the GTP-Schur method, two 
different preconditioners - the local and the Neumann-Neumann preconditioners- are 
considered and are first validated for a simple test case with one fracture. For the 
GTO-Schwarz method, the optimized parameter is used to accelerate the conver¬ 
gence of the associated iterative algorithm. Preliminary numerical experiments show 
that the GTO-Schwarz method converges much faster than the GTP-Schur method 
(with either the preconditioner) in terms of the number of iterations. The Neumann- 
Neumann preconditioner works better than the local preconditioner in the sense that 
its convergence is faster and is only weakly dependent on the mesh size of the dis¬ 
cretizations. The GTO-Schwarz method also has a weak dependence on the mesh size. 
When nonconforming time steps are used, only the local preconditioner preserves the 
accuracy in time: the L 2 error in the fracture of the nonconforming time grid is close 
to that of the conforming fine grid. For the other algorithms, the L 2 error in the 
fracture of the nonconforming time grid is close to that of the conforming coarse grid 
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instead. However, for the GTO-Schwarz method, this weak point when different time 
steps are used is compensated by the fast convergence of the algorithm. 

Appendix A. Proof of Theorem |1.2[ 

We now give the proof of Theorem |1.2| The proof of the theorem is based on the 
Galerkin method, and its main steps will be given after the following lemma, which 
states the main energy estimates. 

Remark A.l. The proof of Lemma \A.2\ is given in the infinite dimensional 
setting but some technical points (those involving u at time t = 0) can only be defined 
using their finite dimensional Galerkin approximations (as was done in detail for 
Dirichlet and Robin boundary conditions in ISSf ). The results presented below have to 
be understood in that sense. 

Lemma A.2. Under assumptions (HI), (H2) and (H4| above, the following a 
priori estimates hold. 


IHl!°“(0,T;M) ^ C 


l U II.L 2 (0,T;E a ) - C 


'tP\\L 2 (0,T-,M) ^ C 
|-® U ll! 2 (0,T;M) ^ C 


^lli 2 (0, T-M) 

^IIl 2 (0,T;M) 

^IIl 2 (0,T;M) 


L\\ 


L 2 (0,T;M) 


I boll 


M 


Iboll 2 


Ml 

Iboll 


M I > 


W 


(A.l) 


w 


where we recall that E a denotes the space E with the norm induced by the bilinear 
form a. 

Proof. As usual we proceed by estimating successively p , u and dtp. 


• First, to derive an estimate for p , we take pit) £ M and u{t) £ E as the test 
functions in © and add the two equations to obtain 

a(u.u) + c(p,p) + ( d t p,p)M = L{p)- 

Using the Cauchy-Schwarz inequality, we see that 


^Ibllif + c (p,p) + NIL < \ (\\ l \\m + \\p\\m) , 


(A.2) 


Now integrating (A.2) over (0, t) for t £ (0, T], we find 

IbWIlL + 2 f c{p,p ) + 2 f ||«HL < IboIlL + II^IIL( 0 .T;M) + / 

JO Jo Jo 


Ibll 


2 

M* 


Then we use the non-negativity of c and apply Gronwall’s lemma to obtain the first 
two estimates in (1 VT1 ) 


Ibll 


L°° (0,T;M) — 


< c 


PoIIm 


m\ 2 L 2 (0,T-,M) 


and 


blll 2 (0,T;E a ) — C 


(I 


PoIIm 


li 2 (0,T;A/) 


(A.3) 


• Next, to estimate dtp , we differentiate the first equation of (1.1) with respect to t 
and take u as a test function. This yields 


a(d t u,u) — b(u,d t p) = 0. 


(A.4) 
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Then taking d t p as a test function in the second equation of (1.1) wee see that 


(d t p, dtP) M + c(p, d t p) + b(u, d t p) = L(d t p). 
Now adding (A.4) and ( |A.5[ ), we obtain 

a(d t u,u) + (d t p, d t p )m + c(p,d t p) = L(d t p), 


(A.5) 


or 


II%>IIm + ^c(p,p) + \j t W u W% a < ^||A||m + \\\ d tP\\ 2 M- (A.6) 

Integrating this inequality over (0, t) for t £ (0, T], we have 

[ II^pIIm + c (p(t)iP(t)) + ll u (0lll o , < \m\ 2 L^(o,T-M) + CcIIpoIIm + ll w (0)||| a j (A.7) 

Jo 

where C c is the constant of continuity of the bilinear form c. There remains to bound 
the term ||M(0)||| a - Toward this end, we use the first equation of <0) with v — u 
and for t = 0: 


a(u(0),u(0)) = b(u(0),p o ). 

The (regularity) assumption that po £ W enables us to write 

ll«( 0 )||s a < C b \\p 0 \\w, 


and, as c(p(t),p(t)) > 0, this gives the third inequality in (A.l) 

• We now derive the last estimate. For this, we take /i = Bu as test function in 
the second equation of 0. 

(d t p, Bu)m + c(p, Bu) + b(u, Bu) = (L, Bu)m 


which we rewrite as 

\\Bu\\ 2 m = (L — d t p, Bu) m - c(p,Bu) <C(\\L\\ 2 m + \\d t p\\ 2 M + \\p\\ 2 M ) +^\\Bu\\ 2 m 


and the fourth inequality then follows by integrating in time and using the previous 
inequalities, which completes the proof of the lemma. □ 


We now give the proof of the theorem. 
Proof. We first prove an estimate for ||u| 


e, which follows easily from the second 


and fourth inequalities in Lemma A.2 and hypothesis (H3’ 


/5|I u IIl 2 (0,T;£) 



WBuWm) ^ C (\\M\ 2 l 2 (o,t- M ) + Ibollvv) ■ 


(A.8) 


Note that this is the only place in the proof where (H3’ 
independant of this hypothesis. 


was used. Lemma A.2 


is 


With the a priori estimates from Lemma A.2 and (A.8), the proof is concluded 
by using Galerkin’s method. □ 
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